Mean field effects on the scattered atoms in condensate collisions 
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We consider the collision of two Bose Einstein condensates at supersonic velocities and focus on the halo of 
scattered atoms. This halo is the most important feature for experiments and is also an excellent testing ground 
for various theoretical approaches. In particular we find that the typical reduced Bogoliubov description, com- 
monly used, is often not accurate in the region of parameters where experiments are performed. Surprisingly, 
besides the halo pair creation terms, one should take into account the evolving mean field of the remaining con- 
densate and on-condensate pair creation. We present examples where the difference is clearly seen, and where 
the reduced description still holds. 



C/2 



o 



> 
m 
m 
in 
in 



X 



I. INTRODUCTION 

When two Bose-Einstein condensates (BECs) collide at 
sufficiently high velocity, pairs of atoms are scattered out 
of the condensates. After many scattering events, a distinct 
halo of atom pairs is formed in momentum space. This was 
observed in many experiments iflUlTIl and analyzed in nu- 
merous theoretical works |[lj, 0, Illlll8l - t34ll . The formation 
of a halo starts spontaneously and is analogous to the gener- 
ation of photon pairs in parametric down conversion. Such 
pho ton pairs were used to observe Bell inequality violation 
ll35ll . and can be applied for quantum cryptography ll36ll or 
quantum teleportation IstIi . In analogy, atoms formed in the 
collision of two BECs have a potential application for preci- 
sion measurements Issll. iri terferometrv Il2l [39l - l42ll . or tests 
of quantum mechanics ll43ll . 

The simplest model that captures the formation of a halo 
in the condensate collision is the "reduced Bogoliubov" 
model (RBM). In this formulation, the two condensates 
counter-propagate at a constant relative velocity and with- 
out change of shape. The wave-packets enter in the Bogoli- 
ubov equation for the field of scattered atoms as a classi- 
cal source iT7y2l|25i|26i|28i[3li[32|l. The RBM can be 
used to calculate various observables, such as the density of 
scattered atoms or their second-order correlation functions, 
which can be directly compared with the experimental data. 
The predictions of the RBM are often in good agreement 
with experiment [31,1 , but there are cases where more com- 
plete models have shown significant departures from the 
RBM flSilll. 

In this paper we compare results of the RBM and the 
complete Bogoliubov equation for a range of parameters. 
In Section we introduce both models and the numerical 
method of solving the equations of motion. Two additional 
simplified formulations, also introduced in Section [III are 
used to investigate the properties of the halo. The depen- 
dence of the halo shape on the model is shown in Sec- 
tion nil Al for two characteristic cases, whereas Sec. HUB I 
explains in detail the dependence of the halo position on the 
degree of simplifications made. We also indicate the range 



of interaction strength for which the RBM is quantitatively 
accurate. 



II. BOGOLIUBOV APPROACH TO THE SCATTERING 
OF ATOMS IN A BEC COLLISION 

A. Scattering in condensate collisions 

Initially the gas is trapped by a harmonic potential. In 
order to start the (half-) collision, a superposition of two 
counter-propagating, mutually coherent, atomic clouds of 
equal density is prepared by shining a Bragg pulse. The 
trap is simultaneously turned off. The two fractions move 
apart with relative speed 2vicc along the z axis, where Wioc 
is the atomic recoil velocity. If the relative velocity is larger 
than the speed of sound at the condensate center, the colli- 
sion is supersonic, superfluidity breaks down and a halo of 
scattered atoms is formed. The physical properties of this 
halo are the main subject of our studies. 

In our approach we model the process of the halo forma- 
tion using a time dependent Bogoliubov method, where the 
field operator is defined as 



^'(x,t) = 0(x,i) +(5(x,t). 



(1) 



Here, 0(x, t) is the condensate wave function governed by 
the GP equation 



ihdt4>{x, t) 



2m 



V2+.g|0(x,t)n0(x,i). (2) 



It is normalized to - the number of atoms in the conden- 
sate, which remains undepleted during the dynamics. The 
coupling constant g relates to the scattering length a through 
g = Anfi^a/m, where m is a mass of an atom. The field op- 
erator (5(x, t) describes the non-condensed atoms (both the 
quantum depletion and the scattered halo) and satisfies the 
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FIG. 1. Cross-sections through the density of non-condensed 
atoms in momentum space nnc(k) — (5^(k)(5(k)) in the kx = 
plane after the end of the collision, for the bose enhanced case of 
a = 27.6, P = 31.6. The labels A, B, C, D correspond to the four 
models introduced in the text. Color scale varies between plots. 
The black circle shows |k| = fco for reference. The red (white) 
contour is at 0.1 (10~^) of the peak condensate density. 



following equation 



ihdt5{x, t) 



2m 

+ ,9 0^(x,t)?t(x,t) 



V2+2,g|0(x,t)|2U(x,t) 



(3) 



From this equation, we can trace back an effective Bogoli- 
ubov Hamiltonian: 



a) 
(4b) 



i?eff = J d3x?t(x) (-|^V2^ ?(X) (4 

+ 2gJ d-^x|0(x)|2?t(x)?(x) 

+ I y" d3x02(x)?^(x)?t(x) + h.c. (4c) 

Line (l4al i stands for kinetic energy, while line ( l4b] l re- 
sults from the interaction of the condensate mean-field 
with the scattered atoms. Finally (l4c] i describes the cre- 
ation/annihilation process of pairs of non-condensed atoms. 

Most salient features of BEC collisions are already found 
for the case of initially spherically symmetric condensates. 
Before the numerical simulation of the collision, we find the 
condensate wave-function as a ground state of the Gross- 
Pitaevskii equation with harmonic trap of frequency u, 



M^) = f-^-V' + imcj^x^ + g\M^)n <^o(x). 
\ 2m 2 / 

(5) 

Next, the trap is turned off and at the same time a set of 
Bragg pulses is applied ifTill . In the center-of-mass frame. 



the initial state is then 

</)(x, 0) = A0o(x) (e*-- + e-*'=«-) , (6) 
where A is the normalization constant and fco = murec/^ is 
the wave-vector associated with the recoil velocity. The ini- 
tial state of the non-condensed part is taken to be a vacuum 



(5(x,0)|0) = 0. 



(7) 



From this initial state the system evolves according to equa- 
tions (|2]i and (O. 



B. Dimensionless parameters 

We point out that in our approximate description (Bo- 
goliubov method), there is a universal scaling such that the 
non-condensed field S{x) is identical in all systems hav- 
ing the same value of gN (or aN, where a is a scatter- 
ing length). Hence the dynamics is described by the length 
scale aN rather than by the scattering length a or the num- 
ber of particles N separately. Other important length scales 
are aho = y/fi/moj - the harmonic oscillator length, and 
1/fco. Alternatively, since there is no trap for t > 0, we 
can use the width of the initial condensate a instead of aho- 
If we apply the Thomas-Fermi approximation ll45ll . a good 
choice of characteristic width a is the Thomas Fermi ra- 
dius Rtf = (15Ng/4:TTmuj^y^^ . In conclusion: there are 
three relevant length scales, hence our system can be char- 
acterized by two dimensionless parameters. We choose: 
/3 = fcoCT and a = aN/a (see 1125;. ^26„ .28., 3 1.1). 

The first parameter /? is independent of interactions and 
has a kinetic character. It can be viewed as the number 
of fringes created by the Bragg pulses on the initial con- 
densate, or as a ratio of the dispersion time scale, nia^ /h 
to the collision timescale, (ma)/ (hko), see for instance 
miHllilll. In reafistic situations ^ > 1. 

The second parameter a is proportional to the ratio of 
the interaction energy per particle gN/ to the kinetic en- 
ergy per particle in the initial condensate f? j {2ma^). It is 
related to a^P''''*^ = moiy/(CT?i^7r'^/^) which has been used 
previously ll25ll26ll28ll3lll . but differs by a numerical factor: 
a = {\/n/4)a^P^^^\ The ratio between a and /3 quantifies 
the strength and nature of scattering. It has been shown that 
Bose enhancement of scattering occurs for a> (3 Ii26i1 . 



C. Simulation of dynamics using the STAB method 

The Bogoliubov equation ([3j can be solved numerically 
using the positive-P representation method lUl, where the 
field operator S is replaced with two complex-number fields 
?/'(x) and tjj{x). The dynamics of the system is governed by 
a pair of stochastic Ito equations. 



J 
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2m 
2m 



V2 + 2g\^i^,t)\^j V(x,t) +g,^(x,i)V(x,t)* + v^0(x, t)e(x, i), (8a) 
V2 + 2g\cl){^,t)\A i>{^,t) + g4>{^,tf^{^,ty + V^0(x, i)e(x, t). (8b) 
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FIG. 2. As Figure^ but for the case of Of = 9.2, (3 — 31.6 where 
bosonic stimulation of scattering is negligible. 



Here ^(x, i) and ^(x,t) are delta-coiTelated, independent, 
real stochastic noise fields with zero mean. The second mo- 
ments are equal to 

(C(x,Oe(x',t')) =0^ and 

(e(x, <)C(x', t')) = (e(x, t)e(x', t')> = - - t'). 

Numerically, ^ and ^ are approximated by real gaussian ran- 
dom variables of variance l/{AtAV) that are independent 
at each point at the computational lattice (of volume AV), 
and at each time step of length At. 

Any physical quantity is obtained by substituting (5^ — )■ 
V'* and S ^ ip and changing from a quantum average of the 
normally-ordered operator to a stochastic average ll46ll . 



The braces (■)st denote the statistical average over re- 
alizations. Observables in k-space follow from the Fourier 
transformation. For example, the one-particle density reads 

Pi(k,k') = <^(k)Xk') +Re(V^(k)Xk')>.t. 

The stochastic equations ^ strictly reproduce the full 
quantum dynamics described by H^ff when the number of 
samples tends to infinity. For finite sample sizes, one ob- 
tains an estimator for the full quantum dynamics, with an 
uncertainty that is calculated by standard methods ll47ll . 
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D. Partially reduced Bogoliubov models 



A fortuitous "side-effect" of the method formulated 
above is the possibility to systematically add or remove 
parts of the effective Hamiltonian (|4]i and GP equation (|2]i to 
understand their impact on the dynamics of the 5 field. This 
includes a numerical simulation of the RBM.To do this, we 
first write the condensate wavefunction as a sum of left- and 
right-moving wavepackets, 

0(x,i) =<^L(x,i)+0fl(x,O, (9) 
where (/)i/^(x, 0) = A^o(x)e^*'^"^ with the GP equation 

— V2 (10a) 
2m 

+5|0L/fl(x,t)P (10b) 

+2.9|0H/i(x,O|' (10c) 

+5 0fl,/i(x,t)*0i/flXx,t) (lOd) 

0L/fl(x,i)- 

The terms proportional to the coupling constant g can be 
interpreted as the self-interaction of the wavepacket ( llObb 
and cross-interaction between different wavepackets ( 11 Oct 
respectively. The remaining processes are in line dlOdt . 

Next, the decomposition (|9]l is put into the Bogoliubov 
Hamiltonian ^ giving 



eff 



d3x?t(x) 



2m 



5(x) 



(11a) 



2g / d^x |0l(x) + ,^fi(x)|'?t(x)?(x) (lib) 



d3x0i(x)0H(x)<5t(x),5t(x) 



h.c. (11c) 



9 



d^X (</)i(x)2 +0^(x)2)?t(x)?t(x) + h.c. (lid) 



The resonant term ( II Icb governs the process where two 
atoms, one from 0/f(x) and one from (/)l(x), collide and 
elastically scatter into the halo localized around the radius 
|k| « /co as in the usual RBM. The line dl Idl l leads to cre- 
ation of two atoms in the 5 field originating from one wave- 
packet, and we call this here "on-condensate pairing". This 
incudes quantum depletion processes. 

Note that the specific form of Equations (fToUTTT i allows 
for removal of chosen terms. This way, one can inspect their 
role in the dynamics of the system. For instance, we can ne- 
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FIG. 3. Cross-sections through the halo of scattered atoms along 
the ky axis after the end of the collision for a — 27.6 (a) and 
a = 9.2 (b). The curves are calculated using the models A (red), 
B (green), C (black) and D (blue). Triple lines indicate the best sta- 
tistical estimate and the ±1(t statistical uncertainties in the mean. 



gleet all terms but the kinetic energy ( 1 11 al l and pair produc- 
tion jllcb in the Bogoliubov Hamiltonian. This simplifica- 
tion we call the Pair Production (PP) dynamics. Within this 
approximation, the Positive P formulation takes the form: 

ihdtip{x,t) = — tp{x,t) + 2g(j)L{y^,t)(j)R{-K,t)'4>{-K,ty 



2m 



(12a) 



+ \/2ihg(t>L{yi, t)(l)R{-x, t) ^(x, t), 
2m 



+ V2«%<^L(x,t)0fl,(x,t)e(x,t). (12b) 

We can also simplify the GP equation (fTTT i neglecting all 
the nonlinear terms and approximating the kinetic energy 
operator As a result of this treatment we obtain the stiff 
movement (SM) of the counter-propagating wave-packets, 

ihdt(l)Li^,t) = (ko ~ 2%^ 0L(x,t) (13a) 



ihdtit'B{x,t) 



2m 
h^kQ 
2m 



ko 



dz 

2z^)0H(x,t). (13b) 



The RBM, described in literature, consists of both 
the PP and SM simplifications. Also, using alternative 
combinations of the approximations described above, we 
can generate four different models: 

A: GP equation ^ + full Bogoliubov (|8j 

B: SM (O recombined to (|9]n- full Bogoliubov dSJ. 

C: SM dnj H- PP ([nil equation we call it RBM 

D: GP equation (reduced as in ( fTOb ) -i- PP ( fT2] i equation. 

In the following section we compare the distributions of 
scattered atoms obtained using the four above models. 



III. RESULTS 
A. The shape of the halo 

We first investigate the dynamics of the system in two 
characteristic cases, varying only a, while /3 = 31.6 is kept 
constant. It has been previously shown, that for Gaussian 
colliding clouds, significant Bose enhancement of scatter- 
ing occurs when a'P'^'^*^//? « 2 or greater, i.e. a > /3. 
Figures [T] and 12] show the cross-section through the halo of 
scattered atoms in the k^ = plane at the end of the colli- 
sion, for the case of appreciable {a = 27.6) and negligible 
(a = 9.2) bosonic enhancement of scattering. Density pro- 
files of the halo in these two cases are shown in Fig.|3] 

When the dynamics of the 5 field is described by the sim- 
plified models C and D, the density of scattered atoms is 
spherically symmetric. On the other hand, this symmetry 
is lost within models A and B, where the halo of atoms 
is weakened near the condensates. This effect can be at- 
tributed primarily to the mean field dl Ibt . with some en- 
hancement resulting from (II Idl i. Notice that for larger in- 
teraction strength - and consequently larger a - these phe- 
nomena are more pronounced. Moreover, models A and B 
predict some non-condensed atoms appearing on top of the 
BECs (in the same location as the BECs). The latter effect 
results from the term ( 11 Idl i only, and it as related to quan- 
tum depletion. 



B. Position of the halo center 

In the momentum space we denote location of the max- 
imum of the halo density by fcmax- This parameter, briefly 
discussed in Ol, varies between models A-D, as can be seen 
in Figs.[T]and|2]and Fig.|3]in more detail. The shift of the po- 
sition of the halo can be explained when the BECs are mod- 
elled with two counter-propagating plane- waves iHllil], and 
using energy conservation argument. The energy of a par- 
ticle released from the condensates depends on the form of 
the GP equation we use to describe it. The full GP equation, 

with n - the mean density of the system, gives + |.gn, 

while the SM ( [13l l gives a purely kinetic value . On the 
other hand the energy needed to place a particle in a non- 
condensed mode with momentum hk is equal to -^^^ + 2gn 

1-2 7 2 

in the case of full Bogoliubov and in the case of the PP 
equation ( fTSI ). Hence, using energy conservation we obtain 



A: 
B: 
C: 
D: 



2m 



2m 



2m ~ 2m 



2m 



h^k^ 
2m 



2m 



h^k„ 
2m 
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FIG. 4. Dependence of the halo peak position fcmax on a/fi^, cal- 
culated with A (red circles), B (green squares), C (black triangles) 
and D (blue x). Error bars indicate the ±1(7 statistical uncertain- 
ties in the mean. 



This simplified model gives quantitative agreement with the 
results presented in Figure |3] 

Notice that in all models (fc/fco)^ is a linear function of 
the parameter gn In order to generalize this depen- 

dence to the case of non-uniform condensates, we replace 
n with half of the maximal density, which in the Thomas- 
Fermi approximation is n„iax = {15/8TT)N/Rj,p. Recall- 
ing the definitions of a and /3, we obtain the shift propor- 
tional to a/ 0^. To verify this conjecture, in Fig|4]we plot 
(fcmax/fco)^ as a function of a/ 0^ for all four models. We 
observe that for growing a//3^ (and thus growing a, as we 
keep constant /3 = 31.6), there is some deviation from the 
linear behaviour for models A and D. What these two mod- 
els have in common is the full evolution of the condensates, 
which might make the picture of plane wave collisions with 
a time-independent density difficult to uphold. 



Finally, we note that for small enough a, the discrepancy 
between the RBM and full Bogoliubov treatment is negligi- 
ble, as expected. The region of agreement (a < (3) matches 
approximately the region where stimulated scattering (Bose 
enhancement) is absent, although investigation with various 
values of (3 would be necessary to confirm or refute a link. 



IV. CONCLUSIONS 

We have demonstrated how different approximations of 
the full Bogoliubov equation dramatically influence the re- 
sults of numerical simulations of BEC collisions. Although 
the simpliest RBM approach predicts a spherical density 
of scattered atoms, more involved models show significant 
discrepancies from this distribution. These differences can 
be understood as a result of the action of mean-field terms 
in the GP and Bogoliubov equations. Also, the position 
of the density maximum changes from model to model. 
This effect, which was investigated when the BECs are ap- 
proximated by the plane-waves, is here interpreted for non- 
uniform condensates. The RBM model is quantitatively 
correct when the interaction strength, as quantified by a, 
is sufficiently small. 
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